
#define K_LOOP_SIZE (SIZE - 1)

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <assert.h>

#define SIZE 4096

double M[4096 * 4096];
;
double L[4096 * 4096];
;
double Momp[4096 * 4096];
;



double t() {
    struct timeval tv;
    assert (gettimeofday(&tv, ((void *)0)) == 0);
    return (double)tv.tv_sec + (double)tv.tv_usec/1000000.0;
}

void lu_llc(int nthreads)
{
    int i;
    int j;
    int k;
    double time_start, time_end;
    for (i = 0; i < SIZE; i++) {
            for (j = 0; j < SIZE; j++) {
                    if (i == j)
                            M[(j * 4096) + i] = 1.0;
                    else
                            M[(j * 4096) + i] = (1.0 + (i * SIZE)) + j;
                    L[(j * 4096) + i] = 0.0;
            }
    }

    ;
    time_start = t();



#pragma llc context name("lu") copy_in(M, L, i, j, k) copy_out(M)
   for (k = 0; k < K_LOOP_SIZE; k++)  {
        #pragma llc for shared(M, L,  k) private(i, j) // schedule(static, 1)
        for (i = k + 1; i < SIZE; i++) {
            L[(k * 4096) + i] = M[(k * 4096) + i] / M[(k * 4096) + k];
                for (j = k + 1; j < SIZE; j++) {
                    M[(j * 4096) + i] = M[(j * 4096) + i] - (L[(k * 4096) + i] * M[(j * 4096) + k]);
                }
        }
    }

    ;
    time_end = t() - time_start;
    printf("** Lu LLC time %g \n", time_end);
}


void lu_openmp(int nthreads)
{
    int i;
    int j;
    int k;
    double time_start, time_end;
    for (i = 0; i < SIZE; i++) {
            for (j = 0; j < SIZE; j++) {
                    if (i == j)
                            Momp[(j * 4096) + i] = 1.0;
                    else
                            Momp[(j * 4096) + i] = (1.0 + (i * SIZE)) + j;
                    L[(j * 4096) + i] = 0.0;
            }
    }

    ;
    time_start = t();



#pragma omp parallel for shared(M, L) private(i, j) //schedule(static,1)
   for (k = 0; k < K_LOOP_SIZE; k++)  {
        for (i = k + 1; i < SIZE; i++) {
            L[(k * 4096) + i] = Momp[(k * 4096) + i] / Momp[(k * 4096) + k];
                for (j = k + 1; j < SIZE; j++) {
                    Momp[(j * 4096) + i] = Momp[(j * 4096) + i] - (L[(k * 4096) + i] * Momp[(j * 4096) + k]);
                }
        }
    }

    ;
    time_end = t() - time_start;
    printf("** Lu OpenMP time %g \n", time_end);
}

int main() {
   lu_llc(1);
   lu_openmp(1);
   if (different())
    printf("*** Result Incorrect \n");
  return 1;
}

int compare_double(double f1, double f2)
{
  //double precision = 1e-15;
  double precision = 1e-5;

  if (((f1 - precision) < f2) &&
     ((f1 + precision) > f2))
  {
   return 1;
  }
  else
  {
    return 0;
  }
}


int different() {
   int i, j;
   for (i = 0; i < SIZE; i++)
    for (j = 0; j < SIZE; j++)
        if (compare_double(Momp[(j * 4096) + i],M[(j * 4096) + i]))
            return 1;

    return 0;
}

